numba を用いた並列計算
numba について
numba は、 JIT(Just In Time)コンパイラ という技術を利用したPythonの拡張ライブラリで、 Anaconda Python をリリースしている Continuum Analytics, Inc が開発しています。 numba には次のような特徴があります。
関数にデコレータをつけるだけでJITコンパイルされて実行できる
Numpyの配列や関数を用いた処理を高速化できる
CPU/GPUを用いた並列計算ができる
numba は次の環境で動作します。
サポートOS:Windows、OSX、Linux
アーキテクチャ:x86、x86_64、ppc64le、armv7l、armv8l(aarch64)
GPU:NvidiaCUDA。 AMDROC
CPython
NumPy1.15以降をサポート
numba のインストール
numba は次のようにインストールします。
code: condaの場合
$ conda install numba
code: pipの場合
$ pip install numba
Numba の依存モジュールは最小限に抑えられていますが、次のパッケージをインストールすることで、追加の機能が提供されます。
scipy:numpy.linalg()関数のコンパイルができるようになる
colorama:バックトレース/エラーメッセージをカラーでの強調表示が有効になる
pyyaml:YAMLフォーマットのファイルを利用してNumbaの構成することができる
icc_rt:高性能のベクトル数学ライブラリIntel SVMLの使用を許可する(x86_64のみ)
numba の利用方法
numba の使い方は驚くほど簡単で、numba モジュールをインポートしておき、 デコレータ@jit を関数の前に記述するだけです。
Numbaは、デコレートされた関数が呼び出されると、実行のためにJITマシンコードにコンパイルして、コードの全部または一部をネイティブのマシンコードで実行できるようにします。
まず、次のようなコードがあるとします。
code: 01_pure_python.py
import numpy as np
x = np.arange(100).reshape(10, 10)
def go_fast(a):
trace = 0.0
for i in range(a.shape0): return a + trace
print(go_fast(x))
これを numba で高速化した場合の利用例です。
code: 02_numba.py
from numba import jit
import numpy as np
x = np.arange(100).reshape(10, 10)
# 最高のパフォーマンスを得るためにnopythonモードを設定
# これは、@njit() でデコレートすることと同じ
@jit(nopython=True)
def go_fast(a): # 関数は、最初に呼び出されたときにマシンコードにコンパイルされる
trace = 0.0
for i in range(a.shape0): # Numba はループ処理に向いている trace += np.tanh(ai, i) # Numba は NumPy の関数をサポート return a + trace # Numba NumPy の broadcasting をサポート
print(go_fast(x))
code: ipython
In 1: %time %run 01_pure_python.py CPU times: user 155 ms, sys: 44.1 ms, total: 199 ms
Wall time: 201 ms
In 2: %time %run 02_numba.py CPU times: user 877 ms, sys: 166 ms, total: 1.04 s
Wall time: 954 ms
In 3: %time print(go_fast(x)) CPU times: user 1.77 ms, sys: 72 µs, total: 1.85 ms
Wall time: 1.85 ms
numba は初回の関数コールでJITコンパイルが実行されます。そのため、初回の呼び出しはPythonだけのコードより遅くなってしまいますが、2回目以降は108倍高速になっていることがわかります。
numba が向かないケース
前述の01_pure_python.py は numba に向いたコードでした。次のコードは、逆に numba ではうまく処理することができません。
code: python
from numba import jit
import pandas as pd
@jit
def use_pandas(a): # numba に向かない関数
df = pd.DataFrame.from_dict(a) # Numba pd.DataFrame はサポートしていない
df += 1
return df.cov()
print(use_pandas(x))
numba は Pandas の DataFrame をうまく処理できないため、use_pandas() をJITコンパイルすることができず、Python インタプリターにファールバックされて実行されます。この結果、Numbaの内部オーバーヘッドのコストが追加されることになります。
code: bash
$ python 03_panda_dataframe.py
/home/iisaka/Class.Parallel/03_03_Numba/03_panda_dataframe.py:6: NumbaWarning:
Compilation is falling back to object mode WITH looplifting enabled because Function "use_pandas" failed type inference due to: non-precise type pyobject
During: typing of argument at /home/iisaka/Class.Parallel/03_03_Numba/03_panda_dataframe.py (8)
File "03_panda_dataframe.py", line 8:
def use_pandas(a): # numba に向かない関数
df = pd.DataFrame.from_dict(a) # Numba pd.DataFrame はサポートしていない
^
@jit
/home/iisaka/.conda/envs/class_parallel/lib/python3.9/site-packages/numba/core/object_mode_passes.py:151: NumbaWarning: Function "use_pandas" was compiled in object mode without forceobj=True.
File "03_panda_dataframe.py", line 7:
@jit
def use_pandas(a): # numba に向かない関数
^
warnings.warn(errors.NumbaWarning(warn_msg,
/home/iisaka/.conda/envs/class_parallel/lib/python3.9/site-packages/numba/core/object_mode_passes.py:161: NumbaDeprecationWarning:
Fall-back from the nopython compilation path to the object mode compilation path has been detected, this is deprecated behaviour.
File "03_panda_dataframe.py", line 7:
@jit
def use_pandas(a): # numba に向かない関数
^
warnings.warn(errors.NumbaDeprecationWarning(msg,
a b
a 1.0 10.0
b 10.0 100.0
Numba がどのように機能するのか
Numbaは、デコレートされた関数のPythonバイトコードを読み取り、これを関数への引数の型に関する情報と組み合わせます。 コードを分析/最適化して、最後にLLVMコンパイラ/ライブラリを使用して、関数をCPU機能に合わせたマシンコードに変換したバージョンを生成します。 このコンパイルされたバージョンは、関数が呼び出されるたびに使用されます。
numba での自動並列化
numba で @jit(parallel=True) のようにデコレートした関数は、自動的に並列化されて他の最適化を処理する Numba 変換パスが有効になります
ユーザー定義の関数での一部の操作、例えば 配列にスカラー値を追加するような、並列セマンティクス(並列化することができる構文)を持つことがあります。 ユーザープログラムには、そうした操作が多数含まれている可能性があり、それぞれの操作を個別に並列化することもできます。しかし、このようなアプローチでは、キャッシュミスが発生する場合があり、パフォーマンスが低下することがよくあります。
numba の自動並列化では、Numbaはユーザープログラムでそのような操作を識別し、隣接する操作を融合して、自動的に並列で実行される1つ以上のカーネルを形成しようとします。 このプロセスは、ユーザープログラムを変更することなく自動的に生成されます。
並列カーネルを作成するために、Numbaのvectorize() や guvectorize() を使って、手動で生成することもできます。
サポートされている最適化
並列セマンティクスを持ち、並列化を試みるすべての配列操作には次のものがあります。
サポートされているすべてのnumba配列操作:
配列式、Numpy配列間、配列とスカラー間、およびNumpy.ufuncs 間の一般的な算術関数が含まれます。これらは、要素ごとまたはポイントごとの配列演算と呼ばれることがよくあります。
単項演算子:+, -, ~
二項演算子:+, -, *, /, //, ?, %, >>, ^, <<, &, **, |
比較演算子:==, !, <, =<, >, =>
nopythonモードでサポートされているNumpy.ufuncs
vectorize()を介してユーザー定義のDUFunc
NumPy縮約関数(Numpy reduction functions)
sum, prod, min, max, argmin, argmax.
NumPyの配列数学関数
mean, var, std.
NumPy配列作成関数
zeros, ones, arange, linspace
いくつかのランダム関数
rand, randn, ranf, random_sample, sample, random, standard_normal, chisquare, weibull, power, geometric, exponential, poisson, rayleigh, normal, uniform, beta , binomamial, f, gamma, lognomal, laplace, laplace, randint, triangular
行列とベクトル、または2つのベクトルの間のNumpyドット関数
オペランドの次元とサイズが一致する場合、上記の操作では多次元配列
次元またはサイズが混在する配列間でのNumpyブロードキャストの完全なセマンティクスはサポートされておらず、選択された次元全体での削減もサポートされていません。
ターゲットがスライスまたはブール配列を使用した配列選択
割り当てられる値がスカラーまたはスライス範囲またはビット配列に互換性があると推測される別の選択である配列割り当て
functoolsのreduce演算子
明示的な並列ループ
parallel=Trueを与えて有効になったコード変換パスのもう1つの機能は、明示的な並列ループのサポートです。 range()の代わりにNumbaのprange()を使用して、ループを並列化できることを指定できます。 ユーザーは、サポートされている削減を除いて、ループに反復間の依存関係がないことを確認する必要があります。
変数がループ本体の以前の値を使用して二項関数/演算子によって更新された場合、削減は自動的に推測されます。 削減の初期値は、+=、-=、*=、/= 演算子に対して自動的に推測されます。 他の関数/演算子の場合、リダクション変数は、prange()ループに入る直前にID値を保持する必要があります。 この方法での削減は、スカラーおよび任意の次元の配列でサポートされています。
次のコードは1次元での例です。
code: 05_prange_numba.py
import numpy as np
from numba import njit, prange
@njit(parallel=True)
def prange_test(A):
s = 0.0
for i in prange(A.shape0): return s
a = np.array(range(10000))
print(prange_test(a))
2次元の場合は次のようになります。
code: 06_prange_2d.py
from numba import njit, prange
import numpy as np
@njit(parallel=True)
def two_d_array_reduction_prod(n):
shp = (13, 17)
result1 = 2 * np.ones(shp, np.int_)
tmp = 2 * np.ones_like(result1)
for i in prange(n):
result1 *= tmp
return result1
これらの prange()の効果を見てみましょう。
code: python
In 2: # %load 05_prange_1d.py ...: import numpy as np
...: from numba import njit, prange
...:
...: @njit(parallel=True)
...: def prange_test1(A):
...: s = 0.0
...: for i in prange(A.shape0): ...: return s
...:
...: def prange_test2(A):
...: s = 0.0
...: for i in prange(A.shape0): ...: return s
...:
...:
...:
...: a = np.array(range(10000))
...: # %time b = prange_test1(a)
...: # %time b = prange_test1(a)
...: # %time b = prange_test1(a)
...: # %time b = prange_test2(a)
...:
In 3: %time b = prange_test1(a) CPU times: user 443 ms, sys: 27.9 ms, total: 471 ms
Wall time: 472 ms
In 4: %time b = prange_test1(a) CPU times: user 359 µs, sys: 30 µs, total: 389 µs
Wall time: 393 µs
In 5: %time b = prange_test1(a) CPU times: user 0 ns, sys: 405 µs, total: 405 µs
Wall time: 417 µs
In 6: %time b = prange_test2(a) CPU times: user 2.14 ms, sys: 0 ns, total: 2.14 ms
Wall time: 2.14 ms
code: python
In 2: # %load 06_prange_2d.py ...: from numba import njit, prange
...: import numpy as np
...:
...: @njit(parallel=True)
...: def two_d_array_reduction_prod1(n):
...: shp = (13, 17)
...: result1 = 2 * np.ones(shp, np.int_)
...: tmp = 2 * np.ones_like(result1)
...:
...: for i in prange(n):
...: result1 *= tmp
...:
...: return result1
...:
...: def two_d_array_reduction_prod2(n):
...: shp = (13, 17)
...: result1 = 2 * np.ones(shp, np.int_)
...: tmp = 2 * np.ones_like(result1)
...:
...: for i in prange(n):
...: result1 *= tmp
...:
...: return result1
...:
...: # %time a = two_d_array_reduction_prod1(10000)
...: # %time a = two_d_array_reduction_prod1(10000)
...: # %time a = two_d_array_reduction_prod1(10000)
...: # %time a = two_d_array_reduction_prod2(10000)
...:
In 3: %time a = two_d_array_reduction_prod1(10000) CPU times: user 2.84 s, sys: 63.8 ms, total: 2.9 s
Wall time: 2.91 s
In 4: %time a = two_d_array_reduction_prod1(10000) CPU times: user 3.92 ms, sys: 0 ns, total: 3.92 ms
Wall time: 4.17 ms
In 5: %time a = two_d_array_reduction_prod1(10000) CPU times: user 2.73 ms, sys: 0 ns, total: 2.73 ms
Wall time: 2.85 ms
In 6: %time a = two_d_array_reduction_prod2(10000) CPU times: user 13.6 ms, sys: 0 ns, total: 13.6 ms
Wall time: 13.8 ms
どちらも初回はJITでコンパイルする時間が必要になるのでオリジナルより遅くなりますが、2回目以降は高速になっています。
prange()と競合状態
スライスまたはインデックスで指定された要素が複数の並列スレッドによって同時に書き込まれる場合は、配列のスライスまたは要素に縮小するときに注意が必要です。 コンパイラがそのようなケースを検出しない場合、競合状態(Race condition)が発生します。
次のコードは、並列forループを実行したときに、競合状態が誤った戻り値をもたらす場合の例です。
code: python
from numba import njit, prange
import numpy as np
@njit(parallel=True)
def prange_wrong_result(x):
y = np.zeros(4)
for i in prange(n):
return y
この場合、競合するかどうかはCPU(コア)数にも依存します。この資料作成時での numba 0.54.1 では、multiprocessing モジュールの cpu_count() の値が numba.config.NUMBA_DEFAULT_NUM_THREADSにも使用されています。
1CPUのプラットフォームで実行すると期待どおりに動作します。
code: python
In 2: # %load 07_bad_loop.py ...: from numba import njit, prange
...: import numpy as np
...:
...: @njit(parallel=True)
...: def prange_wrong_result(x):
...: y = np.zeros(4)
...: for i in prange(n):
...:
...: return y
...:
...: a = np.ones(10000)
...: # prange_wrong_result(a)
...:
In 3: prange_wrong_result(a) In 4: import multiprocessing In 5: multiprocessing.cpu_count() しかし、CPUが複数個あるプラットフォームでは prange()が並列で動作するため、競合が発生してしまいます。
code: python
In 2: # %load 07_bad_loop.py ...: from numba import njit, prange
...: import numpy as np
...:
...: @njit(parallel=True)
...: def prange_wrong_result(x):
...: y = np.zeros(4)
...: for i in prange(n):
...:
...: return y
...:
...: a = np.ones(10000)
...: # prange_wrong_result(a)
...:
...: # import numba
...: # numba.get_num_threads()
...: # numba.set_num_threads(1)
...: # prange_wrong_result(a)
...:
In 3: prange_wrong_result(a) In 5: numba.get_num_threads() In 6: numba.set_num_threads(1) In 7: prange_wrong_result(a) 競合が発生するコードは次の箇所です。
code: pythooon
これは次のコードとほぼ同等になります。
code: python
演算子 += は配列を読み取り、加算して、代入を行うため、y[0] += x[i] は内部では次のように処理されます。
code: python
ここで、prange()により、同じ配列を共有し、共有配列に読み取りと割り当てを行う複数のスレッド/プロセスが生成されます。ループ本体は複数のスレッド/プロセスによって同時に実行されるため、そこで競合状態が発生します。
どのスレッド/プロセスがいつ実行されるかは単純に非決定的であるため、返される配列に間違った値が含まれ、各要素が異なる可能性があります。そのため、1つの要素に競合状態が発生する場合もあれば、まったく発生しない場合もあり、複数の要素に競合状態が発生する場合もあります。
一方、配列全体の累積を実行することは問題ありません。
code: python
from numba import njit, prange
import numpy as np
@njit(parallel=True)
def prange_ok_result_whole_arr(x):
y = np.zeros(4)
for i in prange(n):
return y
code: python
In 2: # %load 08_whole_array.py ...: from numba import njit, prange, get_num_threads
...: import numpy as np
...:
...: @njit(parallel=True)
...: def prange_wrong_result(x):
...: y = np.zeros(4)
...: for i in prange(n):
...:
...: return y
...:
...: a = np.ones(10000)
...: # prange_wrong_result(a)
...: # get_num_threads()
...:
In 3: prange_wrong_result(a) 並列ループの外側でスライス参照を作成する場合と同様に、次のようになります。
code: python
from numba import njit, prange
import numpy as np
@njit(parallel=True)
def prange_ok_result_outer_slice(x):
y = np.zeros(4)
for i in prange(n):
return y
code: python
In 2: # %load 09_outer_slice.py ...: from numba import njit, prange, get_num_threads
...: import numpy as np
...:
...: @njit(parallel=True)
...: def prange_ok_result_outer_slice(x):
...: y = np.zeros(4)
...: for i in prange(n):
...: return y
...:
...: a = np.ones(10000)
...: # prange_ok_result_outer_slice(a)
...: # get_num_threads()
...:
In 3: prange_ok_result_outer_slice(a) numba は、いくつかの縮約(reductions)をサポートして、競合状態が発生させないように実装しています。
それらの1つが y += です。ここで重要なのは、変数のスライス/要素ではなく、変数自体であるということです。その場合、numbaは非常に賢いことをします。スレッド/プロセスごとに変数の初期値をコピーしてから、そのコピーを操作します。並列ループが終了した後、コピーされた値を合計します。2番目の例を取り上げ、2つのプロセスを使用した場合、おおよそ次のようになります。
code: python
y = np.zeros(4)
y_1 = y.copy()
y_2 = y.copy()
for i in nb.prange(n):
if is_process_1:
if is_process_2:
y += y_1
y += y_2
性能比較
ここでは、2つのコードを例に、最適化の方法による性能を比較してみます。
三角恒等式 $ cos(x)^ 2 + sin(x)^ 2の計算
ベクトルの単純な要素ごとの平方根の合計
すべてのパフォーマンス数値はあくまでも目安であり、特に明記されていない限り、
np.arange(1.e7)の入力でIntel i7-4790 CPU(4Core)で実行されたものです。
No Pythonモード vs オブジェクトモード
一般的なパターンは、@jitで関数をデコレートすることです。これは、Numbaが提供する最も柔軟なデコレータです。 @jitは、基本的に2つのコンパイルモードを含みます。
最初に、デコレータされた関数をNo Pythonモードでコンパイルしようとします。これが失敗した場合は、オブジェクトモードを使用して関数のコンパイルを再試行します。
オブジェクトモードでループリフティング(Loop Lifting)を使用すると、パフォーマンスをいくらか向上させることができますが、Pythonモードを使用せずに関数をコンパイルすることは、パフォーマンスを向上させるための鍵です。 Pythonモードのみが使用されず、コンパイルが失敗した場合に例外が発生するようにするには、デコレータ@njitおよび @jit(nopython=True) を使用できます。
ループ
NumPyはベクトル演算の使用に関する強力なイディオムを開発しましたが、Numbaでのループにも完全に適用することができます。 C言語やFortran言語に精通しているユーザーの場合、このスタイルでPythonを作成すると、Numbaで問題なく動作します。
Numbaがバックエンドで利用しているLLVMはC言語系統のコンパイラです。
code: 10_numpy_vs_loops.py
from numba import njit
import numpy as np
@njit
def ident_np_jit(x):
return np.cos(x) ** 2 + np.sin(x) ** 2
def ident_np(x):
return np.cos(x) ** 2 + np.sin(x) ** 2
@njit
def ident_loops_jit(x):
r = np.empty_like(x)
n = len(x)
for i in range(n):
ri = np.cos(xi) ** 2 + np.sin(xi) ** 2 return r
def ident_loops(x):
r = np.empty_like(x)
n = len(x)
for i in range(n):
ri = np.cos(xi) ** 2 + np.sin(xi) ** 2 return r
a = np.ones(10000)
# %time ident_np(a)
# %time ident_np_jit(a)
# %time ident_np_jit(a)
# %time ident_np_jit(a)
# %time ident_loops(a)
# %time ident_loops_jit(a)
# %time ident_loops_jit(a)
# %time ident_loops_jit(a)
この例は、@njitでデコレートした場合、ほぼ同じ速度で実行されます。デコレータがないと、ベクトル化された関数は数桁高速になります。
table: 性能比較
関数名 @njit 実行時間
ident_np No 0.581s
ident_np Yes 0.659s
ident_loops No 25.2s
ident_loops Yes 0.670s
Fastmath
アプリケーションによっては、IEEE754への厳密な準拠はそれほど重要ではありません。 その結果、追加のパフォーマンスを得るという観点から、数値の厳密さを緩和することができます。 Numbaでこの動作を実現する方法は、fastmathキーワード引数を使用することです。
code: 11_fastmath.py
from numba import njit
import numpy as np
# fastmathがない場合、このループは厳密な順序で累積する必要がある
@njit(fastmath=False)
def do_sum(A):
acc = 0.
for x in A:
acc += np.sqrt(x)
return acc
def do_sum_nojit(A):
acc = 0.
for x in A:
acc += np.sqrt(x)
return acc
# fastmathを使用すると、浮動小数点の再関連付けが許可されるため、
# 縮約(reduction)をベクトル化できる
@njit(fastmath=True)
def do_sum_fast(A):
acc = 0.
for x in A:
acc += np.sqrt(x)
return acc
a = np.ones(10000)
# %time do_sum_nojit(a)
# %time do_sum(a)
# %time do_sum(a)
# %time do_sum(a)
# %time do_sum_fast(a)
# %time do_sum_fast(a)
# %time do_sum_fast(a)
table: 性能比較
関数名 実行時間
do_sum 35.2 ms
do_sum_fast 17.8 ms
場合によっては、可能な高速計算最適化のサブセットのみに適用したい場合があります。 これは、一連のLLVM fast-mathフラグをfastmathキーワード引数に与えることで実行できます。
code: python
In 2: # %load 12_fastmath_subset.py ...: from numba import njit
...: import numpy as np
...:
...: def add_assoc(x, y):
...: return (x - y) + y
...:
...: a = njit(fastmath=False)(add_assoc)(0, np.inf)
...: assert np.isnan(a)
...: a = njit(fastmath=True) (add_assoc)(0, np.inf)
...: assert a == 0.0
...: a = njit(fastmath={'reassoc', 'nsz'})(add_assoc)(0, np.inf)
...: assert a == 0.0
...: a = njit(fastmath={'reassoc'}) (add_assoc)(0, np.inf)
...: assert np.isnan(a)
...: a = njit(fastmath={'nsz'}) (add_assoc)(0, np.inf)
...: assert np.isnan(a)
...:
自動並列化
コードにNumbaがサポートしている並列化可能な操作が含まれている場合、Numbaは複数のネイティブスレッドで並列に実行されるバージョンをコンパイルすることができます。このバージョンはGlobal Interpriter Lock(GIL) を発生させずに実行することができます。自動並列化はparallelキーワード引数を追加するだけで有効になります。
code: 13_parallel.py
from numba import njit
import numpy as np
@njit(parallel=True)
def ident_parallel(x):
return np.cos(x) ** 2 + np.sin(x) ** 2
def ident_serial(x):
return np.cos(x) ** 2 + np.sin(x) ** 2
a = np.ones(1000000)
# %time ident_serial(a)
# %time ident_parallel(a)
table: 性能比較
関数名 実行時間
ident_parallel 112 ms
parallel=True を与えた場合のこの関数の実行速度は、NumPyの約5倍、標準の@njitの6倍です。
Numba並列実行では、OpenMP と同様の明示的な並列ループ宣言もサポートされています。 ループを並列で実行する必要があることを示すには、numba.prange()関数を使用する必要があります。この関数は、Pythonのrange()関数のように動作し、parallel=Trueが設定されていない場合は、単にrange()のエイリアスとして機能します。 prange()で誘導されたループは、驚異的な並列計算と縮約(Reduction)を行うことができます。 合計が順不同で累積されても安全であると仮定して、reduce over sumの例を再検討すると、nのループはprange()を使用して並列化できます。 さらに、この場合、各スレッドが部分和を計算するため parallel=True を使用することで、アウトオブオーダー実行 が有効であるという仮定がすでに行われているため、fastmath=Trueキーワード引数を問題なく追加することができます。 code: 14_parallel_with_fastmath.py
from numba import njit, prange
import numpy as np
# 各スレッドは独自の部分和を累積でき、
# クロススレッド縮退が実行されてれる
@njit(parallel=True)
def do_sum_parallel(A):
n = len(A)
acc = 0.
for i in prange(n):
return acc
@njit(parallel=True, fastmath=True)
def do_sum_parallel_fast(A):
n = len(A)
acc = 0.
for i in prange(n):
return acc
def do_sum_serial(A):
n = len(A)
acc = 0.
for i in prange(n):
return acc
a = np.ones(1000000)
# %time do_sum_serial(a)
# %time do_sum_parallel(a)
# %time do_sum_parallel(a)
# %time do_sum_parallel(a)
# %time do_sum_parallel_fast(a)
# %time do_sum_parallel_fast(a)
# %time do_sum_parallel_fast(a)
table: 性能比較
関数名 実行時間
do_sum_parallel 9.81 ms
do_sum_parallel_fast 5.37 ms
Intel SVML
インテル社は、コンパイラー組み込み関数として使用できる最適化された、ベクトル算術関数を計算するための SVML(Short Vector Mathematical Library) を提供しています。 icc_rt モジュールが環境に存在する場合、またはSVMLライブラリが利用可能である場合、Numbaは可能な限りSVML組み込み関数を使用するようにLLVMバックエンドを自動的に構成します。 SVMLは、各組み込み関数の高精度バージョンと低精度バージョンの両方を提供し、使用されるバージョンは、fastmathキーワードを使用して決定されます。 デフォルトでは、最下位桁単位(Unit in the last place)が1 ULP以内の精度を使用しますが、fastmath=Trueに設定されている場合は、組み込み関数の低精度バージョンが使用されます(4 ULP以内)。
まず、condaを使用してSVMLを取得します。
code: bash
$ conda install icc_rt
@njitのオプションのさまざまな組み合わせを使用して、SVMLの有無にかかわらず、上から恒等関数の例ident_npを再実行すると、次のパフォーマンス結果が得られます(入力サイズnp.arange(1.e8))。 参
NumPyだけの場合は、この関数は5.84秒で実行されました。
table: 性能比較
@njit キーワード SVML 実行時間
None No 5.95s
None Yes 2.26s
fastmath=True No 5.97s
fastmath=True Yes 1.8s
parallel=True No 1.36s
parallel=True Yes 0.624s
parallel=True, fastmath=True No 1.32s
parallel=True, fastmath=True Yes 0.576s
SVMLがこの関数のパフォーマンスを大幅に向上させることがわかります。SVMLが存在しない場合のfastmathの影響はゼロです。
線形代数
Numbaは、No Pythonモードで numpy.linalg のほとんどをサポートします。 内部実装は、LAPACKおよびBLASライブラリに依存して数値演算を行い、SciPyから必要な関数のバインディングを取得します。 したがって、Numbaを使用してnumpy.linalg関数で良好なパフォーマンスを実現するには、十分に最適化されたLAPACK / BLASライブラリに対して構築されたSciPyを使用する必要があります。 Anacondaディストリビューションの場合、SciPyは高度に最適化されたIntel MKLに対して構築されているため、Numbaでもこのパフォーマンスの恩恵を受けることができます。
スレッドレイヤー
Numbaスレッドレイヤー(Threading Layer) について説明します。これは、並列ターゲットをCPUの使用して並列実行を実行するための内部的に使用されるライブラリです。
@jitおよび@njitでのparallel=True を与える
@vectorizeおよび@guvectorizeで target='parallel' を与える
コードがthreadingまたはmultiprocessingモジュール(または他の種類の並列処理モジュール)を使用していない場合、Numbaに付属するスレッドレイヤーはデフォルトで適切に機能します。特になにもする必要がありません。
利用できるスレッドレイヤー
使用可能なスレッドレイヤーは3つあり、次のように名前が付けられています。
tbb IntelTBBを利用するスレッドレイヤー
omp OpenMPを利用するに裏打ちされたスレッドレイヤー
workqueue シンプルな組み込みのワークシェアリングタスクスケジューラ
実際には、デフォルトで機能することが保証されているスレッドレイヤーはworkqueueです。 ompレイヤーには、適切なOpenMPランタイムライブラリが存在する必要があります。 同様にtbbレイヤーには、Intel TBBライブラリが必要です。
これらは次のようにインストールすることができます。
code: conda
$ conda install tbb
code: pip
$ pip install tbb
注意
Linux系ディストリビューションでの互換性の問題やその他の移植性の懸念により、
pipコマンドでインストールされる Numbadでは、OpenMPスレッドレイヤーが
無効になっています。
スレッドレイヤーの設定
スレッドレイヤーは、環境変数NUMBA_THREADING_LAYERを介して、またはnumba.config.THREADING_LAYERへの割り当てを介して設定されます。スレッドレイヤーを設定するためのプログラムによるアプローチを使用する場合は、並列ターゲットのNumbaベースのコンパイルが発生する前に論理的に発生する必要があります。スレッドレイヤーを選択するには2つのアプローチがあります。1つはさまざまな形式の並列実行で安全なスレッドレイヤーを選択する方法、もう1つはスレッドレイヤー名(tbbなど)を使用して明示的に選択する方法です。
安全な並列実行のためのスレッドレイヤーの選択
並列実行は、基本的に4つの形式のコアPythonライブラリから派生します。
最初の3つは、他の手段による並列実行を使用するコードにも適用されます。
threadingモジュールからのスレッド
multiprocessingモジュールのspawn()を利用してプロセスを生成。(Windowsデフォルト)
multiprocessingモジュールからfork()を利用してプロセスを生成。(Unixデフォルト)
multiprocessingモジュールからforkserverを使用してプロセスを生成。(Unixのみ)
新しいプロセスが生成され要求に応じてこの新しいプロセスからプロセスが作成される。
これらの形式の並列処理で使用されているライブラリは、指定されたパラダイムの下で安全な動作を示す必要があります。その結果、スレッドレイヤーの選択方法は、特定のパラダイムに対して安全で、クロスプラットフォームおよび環境に耐性のある方法でスレッドレイヤーライブラリを選択する方法を提供するように設計されています。
設定メカニズムに提供できるオプションは次のとおりです。
default:特定の安全保証を提供しない。(デフォルト)
safe:forksafeとthreadsafeの両方。
tbbパッケージ(Intel TBBライブラリ)がインストールされている必要がある
forksafe:forksafeライブラリを提供する
threadsafe:threadsafeライブラリを提供する
選択されたスレッドレイヤーを検出するために、並列実行後に関数numba.threading_layer()を呼び出すことができます。
code: 15_thread_layer.py
from numba import config, njit, threading_layer
import numpy as np
# 並列ターゲットのコンパイルの前にスレッドレイヤーを設定する
config.THREADING_LAYER = 'threadsafe'
@njit(parallel=True)
def foo(a, b):
return a + b
x = np.arange(10.)
y = x.copy()
# 関数のコンパイルが強制され、
# スレッドレイヤーが選択されてから並行実行される
foo(x, y)
# 選択したスレッドレイヤーを表示
print("Threading layer chosen: %s" % threading_layer())
TBBがインストールされていないLinuxマシンの場合は、次の実行結果となります。
code: bash
Threading layer chosen: omp
Linux系プラットフォームで利用できるGNU OpenMPはスレッドセーフです。
名前付きスレッドレイヤーの選択
上級ユーザーは、ユースケースに合わせて特定のスレッドレイヤーを選択することをお勧めします。これは、スレッドレイヤー名を設定メカニズムに直接指定することによって行われます。 オプションと要件は次のとおりです。
table: 設定
スレッドレイヤー名 プラットフォーム 必要項目
tbb All tbb パッケージ
Linux GNU OpenMP ライブラリ
omp Windows MS OpenMP ライブラリ
OSX intel-openmp パッケージ
workqueue All None
スレッドレイヤーが正しくロードされない場合、Numbaはこれを検出し、問題を解決する方法についてのヒントを提供します。 Numba診断コマンド numba -s には、現在の環境でのスレッドレイヤーの可用性を__Threading LayerInformation__としてレポートします。
スレッドレイヤーの追記
スレッドレイヤーは、CPython内部およびシステムレベルライブラリとかなり複雑な相互作用を持っています。そのため、次のような注意すべきいくつかの追加事項があります。
Intel TBBライブラリをインストールすると、スレッドレイヤーの選択プロセスで利用できるオプションが大幅に広がります。
Linuxでは、GNU OpenMPランタイムライブラリ(libgomp)がフォークセーフではないため、ompスレッドレイヤーはフォークセーフではありません。 ompスレッドレイヤーを使用しているプログラムでフォークが発生した場合、フォークされた子プロセスを正常に終了し、エラーメッセージを標準エラー出力に出力する検出メカニズムが存在します。
fork()システムコールが使用可能なシステムで、TBBでバックアップされたスレッドレイヤーが使用されていて、TBBを起動したスレッド以外のスレッド(通常はメインスレッド)からフォーク呼び出しが行われると、動作が未定義になり、標準エラー出力に警告が表示されます。 spawn()は基本的にfork()の後にexec()が続くため、メインスレッド以外からspawn()しても安全ですが、これは単なるfork()呼び出しと区別できないため、警告メッセージが引き続き表示されます。
OSXでは、OpenMPベースのスレッドレイヤーを有効にするためにintel-openmpパッケージが必要です。
スレッド数の設定
numbaが使用するスレッドの数は、numba.config.NUMBA_DEFAULT_NUM_THREADSを参照して知ることができ、これはシステムの使用可能なCPUコア数に基づいています。この値は、環境変数NUMBA_NUM_THREADSで上書きすることができます。
numbaが起動するスレッドの総数は、変数numba.config.NUMBA_NUM_THREADSにあります。
いくつかのユースケースでは、スレッド数をより低い値に設定して、numbaをより高いレベルの並列処理で使用できるようにすることが望ましい場合があります。
スレッド数は、numba.set_num_threads()を使用して実行時に動的に設定できます。 ただし、set_num_threads()では、スレッド数をNUMBA_NUM_THREADSよりも小さい値しか設定できないことに注意してください。 Numbaは常にnumba.config.NUMBA_NUM_THREADSの数のスレッドを起動しますが、set_num_threads() により、未使用のスレッドがマスクされ、計算に使用されなくなります。
numbaが使用している現在のスレッド数には、numba.get_num_threads() を使用してアクセスできます。どちらの関数も、JIT関数の内部で機能します。
スレッド数の制限例
この例では、実行しているマシンに8つのコアがあるとします。つまり、numba.config.NUMBA_NUM_THREADSは8。 @njit(parallel=True) を使用してコードを実行しするとき、4つの異なるプロセスでコードを同時に実行したいとします。デフォルトのスレッド数では、各Pythonプロセスは8スレッドを実行するため、合計で4 * 8 = 32スレッドになります。これは、システムのもつ8コアより大きくなってしまいます。この場合、むしろ各プロセスを2スレッドに制限して、合計が4 * 2 = 8になるようにして、物理コア数と一致するようにします。
これを行うには2つの方法があります。ひとつは、環境変数NUMBA_NUM_THREADSを2に設定することです。
code: bash
# この実行だけに有効
$ env NUMBA_NUM_THREADS=2 python ourcode.py
code: bash
$ export NUMBA_NUM_THREADS=2
# これ以後、シェルが終了するか設定が変更されるまで有効
$ python ourcode.py
このアプローチの利点は、コードを変更せずにプロセスの外部から実行できることです。
ただし、この方法には2つの欠点があります。
環境変数NUMBA_NUM_THREADSは、Numbaをインポートする前、理想的にはPythonを起動する前に設定する必要があります。 Numbaがインポートされるとすぐに、環境変数が読み取られ、Numbaが起動するスレッドの数に応じて、そのスレッドの数がロックされます。起動後にプロセスで使用されるスレッドの数を増やすことができません。
環境変数NUMBA_NUM_THREADSは、プロセスに対してのスレッドの最大数を設定します。 numba.config.NUMBA_NUM_THREADSより大きい値でset_num_threads()を呼び出すと、エラーが発生します。
別のアプローチは、コード中で numba.set_num_threads()関数を使用することです。
code: python
from numba import njit, set_num_threads
@njit(parallel=True)
def func():
...
set_num_threads(2)
func()
並列コードを実行する前にset_num_threads(2) を呼び出すと、並列コードが2つのスレッドでのみ実行されるという点で、環境変数NUMBA_NUM_THREADS=2としてプロセスを実行するのと同じ効果があります。ただし、後でset_num_threads(8)を呼び出して、スレッドの数をデフォルトのサイズに戻すことができます。また、Numbaがインポートされる前に設定されているどうかを心配する必要はありません。並列関数を実行する前に呼び出す必要があるだけです。
APIリファレンス
numba.config.NUMBA_NUM_THREADS
numbaによって起動されたスレッドの総数(最大)。
デフォルトはnumba.config.NUMBA_DEFAULT_NUM_THREADSですが、NUMBA_NUM_THREADS環境変数でオーバーライドできます。
numba.config.NUMBA_DEFAULT_NUM_THREADS
システム上のCPUコアの数(multiprocessing.cpu_count()によって決定されます)。 NUMBA_NUM_THREADS環境変数が設定されていない限り、これはnumba.config.NUMBA_NUM_THREADSのデフォルト値です。
numba.set_num_threads(n)
並列実行に使用するスレッド数を設定します。
デフォルトでは、すべてのnumba.config.NUMBA_NUM_THREADSスレッドが使用されます。
この機能は、使用されていないスレッドをマスクすることによって機能します。したがって、スレッド数nは、起動されるスレッドの総数であるNUMBA_NUM_THREADS以下である必要があります。
この関数は、JIT関数内で使用できます。
n:スレッドの数。 1からNUMBA_NUM_THREADSの間でなければなりません。
numba.get_num_threads()
並列実行に使用されるスレッドの数を取得します。
デフォルトでは(set_num_threads()が呼び出されない場合)、すべてのnumba.config.NUMBA_NUM_THREADSスレッドが使用されます。
この数は、起動されるスレッドの総数numba.config.NUMBA_NUM_THREADS以下です。
この関数は、JIT関数内で使用できます。
戻り値:スレッドの数
自動並列化の診断レポート
jit()の並列オプションは、デコレートされた関数を自動的に並列化する際に実行される変換に関する診断情報を生成することができます。 この情報には2つの方法でアクセスできます。
環境変数NUMBA_PARALLEL_DIAGNOSTICSを設定する方法
parallel_diagnostics()を呼び出す方法
どちらのメソッドも同じ情報を提供し、標準出力へ出力します。 診断情報の冗長レベルは、1から4までの値の整数引数によって制御されます。1は最も冗長性が低く、4は最も冗長性が高くなります。
code: 20_parallel_diagnostics.py
from numba import njit, prange
import numpy as np
@njit(parallel=True)
def test(x):
a = np.sin(x)
b = np.cos(a * a)
acc = 0
for i in prange(n - 2):
for j in prange(n - 1):
return acc
a = np.arange(10)
test(a)
test.parallel_diagnostics(level=4)
code: bash
================================================================================
======= Parallel Accelerator Optimizing: Function test, example.py (4) =======
================================================================================
Parallel loop listing for Function test, example.py (4)
--------------------------------------|loop #ID @njit(parallel=True) |
def test(x): |
a = np.sin(x)---------------------| #0 b = np.cos(a * a)-----------------| #1 acc = 0 |
for i in prange(n - 2):-----------| #3 for j in prange(n - 1):-------| #2 return acc |
--------------------------------- Fusing loops ---------------------------------
Attempting fusion of parallel loops (combines loops with similar properties)...
Trying to fuse loops #0 and #1: - fusion succeeded: parallel for-loop #1 is fused into for-loop #0. Trying to fuse loops #0 and #3: - fusion failed: loop dimension mismatched in axis 0. slice(0, x_size0.1, 1)
!= slice(0, $40.4, 1)
----------------------------- Before Optimization ------------------------------
Parallel region 0:
+--0 (parallel)
+--1 (parallel)
Parallel region 1:
+--3 (parallel)
+--2 (parallel)
--------------------------------------------------------------------------------
------------------------------ After Optimization ------------------------------
Parallel region 0:
+--0 (parallel, fused with loop(s): 1)
Parallel region 1:
+--3 (parallel)
+--2 (serial)
Parallel region 0 (loop #0) had 1 loop(s) fused. Parallel region 1 (loop #3) had 0 loop(s) fused and 1 loop(s) serialized as part of the larger parallel loop (#3).
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
---------------------------Loop invariant code motion---------------------------
Instruction hoisting:
Failed to hoist the following:
dependency: $arg_out_var.10 = getitem(value=x, index=$parfor__index_5.99)
dependency: $0.6.11 = getattr(value=$0.5, attr=sin)
dependency: $arg_out_var.17 = $expr_out_var.9 * $expr_out_var.9
dependency: $0.10.20 = getattr(value=$0.9, attr=cos)
Has the following hoisted:
$const58.3 = const(int, 1)
$58.4 = _n_23 - $const58.3
--------------------------------------------------------------------------------
code: python
@njit(parallel=True)
def test(n):
for i in prange(n):
temp = np.zeros((50, 50)) # np.zeros() で一時配列を割り当てる
for j in range(50):
内部的には、これはおおよそ次のように変換されます。
code: python
def test(n):
for i in prange(n):
temp = np.empty((50, 50)) # p.zeros() で一時配列を割り当てる
for j in range(50):
それから、次のように引き上げられます。
code: python
@njit(parallel=True)
def test(n):
# np.empty()は同値で満たされているため、割り当てはループ不変条件として引き上げられます
temp = np.empty((50, 50))
for i in prange(n):
temp: = 0 # 割り当ては副作用であるため、これは残る for j in range(50):
np.zeros()の初期化が割り当てとアサインに分割され、初期化が iのループから引き上げられ、割り当てが1回だけ発生するため、より効率的なコードが生成されることになります。
並列診断レポート
並列診断レポートでは次の項目についてレポートを出力します。
コード注釈
並列セマンティクスが識別および列挙されたループを備えた装飾関数のソースコードが含まれています。 ソースコードの右側にある loop #IDカラムには、識別された並列ループが出力されます。
次の例では、#0 は np.si()n、#1はnp.cos()、#2 と #3 は prange() です。
code: bash
Parallel loop listing for Function test, example.py (4)
--------------------------------------|loop #ID @njit(parallel=True) |
def test(x): |
a = np.sin(x)---------------------| #0 b = np.cos(a * a)-----------------| #1 acc = 0 |
for i in prange(n - 2):-----------| #3 for j in prange(n - 1):-------| #2 return acc |
ループIDは、検出された順序で列挙されますが、ソースに存在する順序と必ずしも同じではないことに注意してください。 さらに、並列変換はループIDのインデックス付けに静的カウンターを使用することにも注意してください。 結果として、ユーザーには見えない内部最適化/変換が行われるために同じカウンターを使用するため、ループIDインデックスが0から開始しない可能性があります。
融合ループ(Fusing loops)
検出されたループを融合する試みについて、成功したループと失敗したループについてレポートします。 融合に失敗した場合、理由が示されます(他のデータへの依存など)。
code: bash
--------------------------------- Fusing loops ---------------------------------
Attempting fusion of parallel loops (combines loops with similar properties)...
Trying to fuse loops #0 and #1: - fusion succeeded: parallel for-loop #1 is fused into for-loop #0. Trying to fuse loops #0 and #3: - fusion failed: loop dimension mismatched in axis 0. slice(0, x_size0.1, 1)
!= slice(0, $40.4, 1)
この場合、ループ#0と#1の融合が試みられて、成功したことがわかります(どちらもxの同じ次元に基づいています)。 ループ#0と#1の融合が成功した後、ループ#0(融合されたループ#1を含んでいる)とループ#3の間で融合が試みられて、ループの寸法が一致しないため、この融合は失敗したことがわかります。
ループ#0はサイズx.shapeであるのに対し、ループ#3はサイズx.shape[0] - 2です。
最適化前
最適化が行われる前のコード内の並列領域の構造を示しますが、最終的な並列領域に関連付けられたループがあります(これは、最適化の前後の出力を直接比較できるようにするためです)。 融合できないループがある場合、複数の並列領域が存在する可能性があります。この場合、各領域内のコードは並列に実行されますが、各並列領域は順次実行されます。
code: bash
Parallel region 0:
+--0 (parallel)
+--1 (parallel)
Parallel region 1:
+--3 (parallel)
+--2 (parallel)
融合ループのレポートでわかるように、コードには必然的に2つの並列領域があります。 最初のループにはループ#0と#1が含まれ、2番目のループには#3と#2が含まれます。最適化がまだ行われていないため、すべてのループは並列としてマークされます。
最適化後
最適化が行われた後のコード内の並列領域の構造をレポートします。 この場合も、並列領域は対応するループとともに列挙されますが、今回は融合またはシリアル化されたループが記録され、要約が表示されます。
code: bash
Parallel region 0:
+--0 (parallel, fused with loop(s): 1)
Parallel region 1:
+--3 (parallel)
+--2 (serial)
Parallel region 0 (loop #0) had 1 loop(s) fused. Parallel region 1 (loop #3) had 0 loop(s) fused and 1 loop(s) serialized as part of the larger parallel loop (#3).
並列領域0(Parallel region 0) にはループ#0が含まれていて、融合ループでレポートされているように、ループ#1はループ#0に融合されていることに注意してください。
並列領域1(Parallel region 1) にはループ#3が含まれており、ループ#2(内側のprange())はループ#3の本体で実行するためにシリアル化されていることにも注意してください。
ループの引き上げ
最適化が行われた後の各ループについて示します。
引き上げに失敗した指示と失敗の理由(依存性/不純)
引き上げられたオペレーション
発生した可能性のある割り当ての引き上げ
code: bash
Instruction hoisting:
Failed to hoist the following:
dependency: $arg_out_var.10 = getitem(value=x, index=$parfor__index_5.99)
dependency: $0.6.11 = getattr(value=$0.5, attr=sin)
dependency: $arg_out_var.17 = $expr_out_var.9 * $expr_out_var.9
dependency: $0.10.20 = getattr(value=$0.9, attr=cos)
Has the following hoisted:
$const58.3 = const(int, 1)
$58.4 = _n_23 - $const58.3
まず注意することは、ここで表示される情報は変換される関数のNumba中間表現(IR:Intermediate Representation)を参照しているため、上級ユーザー向けであるということです。 このコードの場合、式 a * a は、IRの式 $arg_out_var.17 = $expr_out_var.9 * $expr_out_var.9 に部分的に変換されます。これはループ不変条件ではないため、明らかにループ#0から引き上げることはできません。 ループ#3では、式 $const58.3 = const(int, 1) は b[j + 1] から取得されますが、数値1は明らかに定数であるため、ループから引き上げることができます。
コマンドラインインターフェイス
NumbaはPythonパッケージです。通常、Pythonからnumbaをインポートし、Pythonアプリケーションプログラミングインターフェイス(API)を使用します。ただし、Numbaにはコマンドラインインターフェイス(CLI)も付属しています。つまり、Numbaのインストール時にインストールされるツールnumbaコマンドです。
現在、CLIの唯一の目的は、システムとインストールに関する情報をすばやく表示したり、Numbaを使用してPythonスクリプトのデバッグ情報をすばやく取得したりできるようにすることです。
使用法
ターミナルからnumbaコマンドを使用するには、以下で説明するように、numbaに続けて、--helpや-sなどのオプションと引数を使用します。
PATHが適切に構成されていないと、numbaと入力すると、「コマンドが見つかりません」というエラーが発生する場合があります。その場合、次のようにPython から使用することもできます。
code: bash
$ python -m numba
numbaとpython-m numba の2つは同じものです。
IPythonまたはJupyterからnumbaコマンドを使用するには、!numba を使用します。つまり、コマンドの前にエクスクラメーションマーク(!)を付けます。これは、シェルコマンドを実行するための一般的なIPython / Jupyter機能であり、通常のPythonターミナルでは使用できません。
ヘルプメッセージ
利用可能なすべてのオプションを表示するには次のようにコマンドを入力します。
code: bash
$ numba --help
positional arguments:
filename Python source filename
optional arguments:
-h, --help show this help message and exit
--annotate Annotate source
--dump-llvm Print generated llvm assembly
--dump-optimized Dump the optimized llvm assembly
--dump-assembly Dump the LLVM generated assembly
--annotate-html ANNOTATE_HTML
Output source annotation as html
-s, --sysinfo Output system information for bug reporting
インストールの確認
PythonのREPLプロンプトからNumbaをインポートできるはずです。
code: bash
numba --sysinfo(または略してnumba -s)コマンドを実行して、システム機能に関する情報を報告することもできます。
code: bash
System info:
--------------------------------------------------------------------------------
__Time Stamp__
2018-08-28 15:46:24.631054
__Hardware Information__
Machine : x86_64
CPU Name : haswell
CPU Features :
aes avx avx2 bmi bmi2 cmov cx16 f16c fma fsgsbase lzcnt mmx movbe pclmul popcnt
rdrnd sse sse2 sse3 sse4.1 sse4.2 ssse3 xsave xsaveopt
__OS Information__
Platform : Darwin-17.6.0-x86_64-i386-64bit
Release : 17.6.0
System Name : Darwin
Version : Darwin Kernel Version 17.6.0: Tue May 8 15:22:16 PDT 2018; root:xnu-4570.61.1~1/RELEASE_X86_64
OS specific info : 10.13.5 x86_64
__Python Information__
Python Compiler : GCC 4.2.1 Compatible Clang 4.0.1 (tags/RELEASE_401/final)
Python Implementation : CPython
Python Version : 2.7.15
Python Locale : en_US UTF-8
__LLVM information__
LLVM version : 6.0.0
__CUDA Information__
Found 1 CUDA devices
compute capability: 3.0
pci device id: 0
pci bus id: 1
まとめ
numba はコード全体を並列化するものではありませんが、特定のNumPy を使った数値計算のループを簡単に並列化できることは注目すべきことです。
参考